824  HDhPTIVE  FINITE  ELEMENT  METHODS  FOR  PhR'HBOLIC  PARTIAL  1/1 
DIFFERENTIAL  EQUflTI.  .  (U>  RENSSELAER  POLYTECHNIC  INST 
TROY  NY  DEPT  OF  MATHEMATICAL  SCIE. . 

UNCLASSIFIED  J  E  FLAHERTY  ET  RL.  MAY  83  RFOSR-TR-83-06S9  F/G  12/1  NL 


i*... 

y  RESOLUTION  TEST  CHART 

BUREAU  OF  STANDARDS- 1963'A 


UNCIJaSSIFiED 

SECUWITY  CLASSiriCATION  OF  THIS  PAGE  fnhmt  DM0  Enffd} 

REPORT  DOCUMENTATION  PAGE 

'XP0®fl[“‘TK.  8  3-  06  89  f“" 


FINITE  ELEMENT  METHODS  FOB  PARABOLIC 
PARTIAL  DIFFERENTIAL  EQUATIONS 


READ  INSTRUCTIONS 

before  completing  form 


1.  RCCIPIENT'S  catalog  NUMBER 


%.  TYPE  OF  ntPOHT  6  PERtOO  COVEnro 


lotorim 


6.  PERFORMING  ORC.  REPORT  NUMBER 


7.  AUTHOKfsJ 

Joseph  E.  Flaherty,  J.  Michael  Coyle,  Raymond 
Ludwig,  Stephen  F.  Davis 


S.  CONTRACT  OR  grant  numBERCcJ 


AFOSR  80-0192 


9-  PERFORMING  organization  NAME  AND  ADDRESS 

Rensselaer  Polytechnic  Institute 
Department  of  Mathematical  Sciences 
Troy,  N.Y.  12181 


controlling  OFFICE  NAME  AND  ADDRESS 

Air  Force  Office  of  Scientific  Research  (NM) 
Bolling  AFB,  Washington,  D.C.  20332 


10.  PROGRAM  element,  project,  TASK 
AREA  A  WORK  UNIT  NUMOERS 

61102F 

9749-03 


U.  REPORT  DATE 

May  1983 


<9.  NUMBER  OF  PAGES 
21 


4.  MONITORING  AGENCY  NAME  A  AOORESSftf  dl/Uranl  Imn  CatitrolUng  Olllet)  IS.  SECURITY  CLASS,  (ot  thim  report; 

UNCLASSIFIED 


iSa.  declassificationTdowncradinc 
schedule 


is.  distribution  statement  (ot  Ulto  Koporl) 


Approved  for  public  release;  distribution  unlimited. 


17.  distribution  statement  (ot  U)a  aSalract  ontotod  In  Block  70,  It  ditloront  tnm  Roport) 


IS.  SUPPLEMENTARY  NOTES 

to  appear  in  the  Proceedings  of  the  ARO  Workshop  on  Adaptive  Methods  for 
Partial  Differential  Equations,  held  14-16  February  1983,  University  of 
Maryland 


If.  KEY  WORDS  (Continum  on  roi^oroo  midm  i/nocoo«ory  ond  IdonrI/x  numbor^ 

Adaptive  methods,  parabolic  partial  differential  equations,  finite  element 
methods 


20k  ABSTRACT  fCooiSouo  «o  rovoroo  oMi  1/  nocooooc^  bj  block  numbor) 

We  discuss  a  finite  element  method  for  solving  initial-boundary  value 
problems  for  vector  systems  of  partial  differential  equations  in  one  space 
dimension  and  time.  The  method  automatically  adjusts  the  computational  mesh 
as  the  solution  evolves  in  time  so  as  to  approximately  minimize  the  local 
discretization  error.  We  arc  thus  able  to  calculate  accurate  solutions  with 
fewer  elements  than  would  be  necessary  with  a  uniform  mesh. 

(continued  on  back) 


00  I  iAM'n  1473  EDITION  or  t  NOV  AS  IS  OBSOLETE 


UNCLASSIFIED 

SECURITY  classification  OF  THIS  PAGE  (When  Doto  Entrrod) 


"Adaptive  Finite  Klcmont  Metliods 

* 

for  Parabolic  Partial  Differential  Equations"  by  Joseph  E.  Klalierty  ,  J.  Michael 
*  *  *• 

Coyle  ,  Raymond  Ludwig  ,  and  Steplien  F.  Davis 


AFOSR-TR.  8  3-  06  89 


(L^ 


7 '/ 


Abstract.  »»lA?""discuss  a  finite  element  method  for  solvidq  initial¬ 
boundary  value  problems  for  vector  systems  of  partial  differential 
equations  in  one  space  dimension  and  time.  The  method  automatically 
adjusts  the  computational  mesh  as  the  solution  evolves  in/ time  so  as  to 
approximately  minimize  the  local  discretization  error.  >/e  are  thus 
able  to  calculate  accurate  solutions  with  fewer  elements  than  would  he 
necessary  with  a  uniform  mesh. 

72  is  i 

-Owr  overall  method  contalr^  two  distinct  steps:  a  solution  step 
and  a  mesh  selection  step,  ye  solve  the  partial  differential  equa¬ 
tions  using  a  finite  eleraent-Galerkiu  method  on  trapezoidal  space-time- 
elements  with  either  piecewise  linear  or  cubic  Hermite  polynomial 
approximations.  A  variety  of  mesh  selection  strategies  are  discussed 
and  analyzed.  Results  are  presented  for  several  computational 
examples. 

\ 


1.  Introduction .  We  consider  adaptive  finite  element  procedures 
for  finding  numerical  solutions  of  M-dimenstona 1  vector  systems  of 
partial  differential  equations  that  have  the  form 


(1.1)  Lu  :=  Ut  +  f(x,t,u,Ux)  -  (D(x, t, u )Ux)x  =  0,  a  <  X  <  b,  t  >  0, 
subject  to  the  Initial  and  linear  separated  boundary  conditions 


(1.2a)  u(x,0)  =  uq(x),  a  <  x  <  b. 
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(1.2b)  Biu(a,t)  ;=  A,^(t)u(a,t)  +  A,  2(  t  )Ujj(a,  t)  =b^(t), 

(1.2c)  B2u(b,t)  ;=  A2i(t)u(b,t)  +  A22(t)u3j(b,  t)  =b2(t),  t  >  0. 

There  are  )t^  initial  boundary  conditions  (1.2a)  and  ^2  terminal  bound¬ 
ary  conditions  (1.2b).  We  are  primarily  concerned  with  parabolic 
problems  where  D  is  positive  definite  and  =  k2  =  M;  however,  we  do 
not  restrict  ourselves  to  this  case,  but  instead  we  assume  that 
conditions  are  specified  so  that  equations  (1.1)  and  (1.2)  have  an 
isolated  solution. 

Problems  having  the  above  form  arise  in  many  applications  and  our 
ultimate  goal  is  to  create  reliable  and  robust  software  that  will 
solve  a  wide  class  of  them  without  requiring  users  to  supply  numerical 
data  such  as  temporal  and  spatial  step  sizes.  Thus,  we  envision  a 
computer  code  that  will  automatically  discretize  and  solve  (1.1), 

(1.2)  on  a  nonuniform  computational  net  and  attempt  to  meet  a  pre¬ 
scribed  error  tolerance. 

The  key  decisions  that  must  be  made  in  developing  such  a  code  are 
selecting  (1)  spatial  and  temporal  numerical  integration  methods, 

(ii)  error  Indicator  and  estimator  procedures,  and  (iil)  adaptive 
static  and  dynamic  mesh  allocation  and  distribution  techniques. 

We  discretize  and  solve  (1.1),  (1.2)  using  a  finite  element- 
Galerkin  method  on  trapezoidal  space-time  elements*  A  detailed  dis¬ 
cussion  and  analysis  of  this  approach  was  given  in  Davis  [6]  and 
Davis  and  Flaherty  (7],  and  herein  we  shall  only  repeat  (cf.  Section  2) 
those  features  that  are  necessary  to  the  continuity  of  this  paper.  Our 
technique  is  similar  to  that  of  Jamet  and  Bonnerot  [11],  and  we  chose 
it  because  it  is  generally  easier  to  generate  high  order  approximations 
to  partial  differential  equations  on  a  nonuniform  mesh  by  using  finite 
element  methods  than  by  using  finite  difference  methods. 

We  combine  the  error  indication  and  static  and  dynamic  mesh  adapta¬ 
tion  steps  by  moving  a  fixed  number  of  finite  elements  so  that  they 
approximately  minimize  the  discretization  error  per  time  step.  This 
task  is  known  (cf.  de  Boor  [8]  or  Pereyra  and  Sewell  [15])  to  be 
asymptotically  equivalent  to  selecting  a  mesh  that  equidistributes  the 
error,  l.e.,  a  mesh  where  the  error  is  equal  on  every  element. 

Our  approach  is  somewhat  similar  to  the  work  of  Miller  et  al. 
[9,13,14],  except  that  they  couple  the  finite  element  solution  and  mesh 
adaptation  steps  whereas  we  consider  them  as  distinct  phases. 

In  contrast  to  these  two  methods,  Bieterman  and  Babuska  [2,3]  use  a 
finite  element  method  of  lines.  In  this  approach  the  partial  differen¬ 
tial  equations  are  first  discretized  in  space  using  a  Galerkln  method. 
This  yields  a  system  of  ordinary  differential  equations  in  time  which 
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may  be  solved  using  one  of  the  many  available  ordinary  differential 
equations  codes.  They  also  add  or  delete  elements  in  regions  where  the 
spatial  discretization  error  is  estimated  to  be  too  large  or  too  small. 
Their  method  can  potentially  solve  partial  differential  equations  to  a 
prescribed  level  of  accuracy,  whereas  this  is  generally  not  possible 
with  moving  mesh  methods  that  use  a  fixed  number  of  elements.  On  the 
other  hand,  moving  mesh  methods  can  follow  evolving  nonuniformities  and 
very  effectively  reduce  dispersive  errors  (cf.  Hedstroro  and  Rodrigue 
[lOj).  Quite  clearly,  a  code  that  had  both  mesh  moving  capabilities 
and  the  ability  to  add  and  delete  elements  would  be  an  ideal  software 
tool  for  solving  partial  differential  equations.  We  have  been  experi¬ 
menting  with  adding  elements  as  the  temporal  integration  proceeds  In 
example  3  of  Section  4. 

In  Section  3  of  this  paper  we  discuss  several  mesh  equldistribution 
algorithms  and  explore  their  properties.  We,  unfortunately,  show  that 
many  algorithms  that  are  based  on  integrating  ordinary  differential 
equations  for  the  mesh  velocities  are  unstable  for  dissipative  partial 
differential  equations.  In  Section  4,  we  apply  our  methods  to  several 
examples  and  in  Section  5  we  discuss  our  results  and  suggest  some 
extensions  and  improvements. 

2.  Finite  Element  Solution  Procedure.  We  discretize  problem  (1.1)- 
(1.2)  on  the  strip 

(2.1)  Sn  :=  {(x,t)  |  a<x<b,  tn<t<  , 

using  a  finite  element-Galerkin  procedure.  Hence,  we  approximate 
u(x,t)  on  Sfj  by  U(x,t)  €  select  "test"  functions  V(x,  t)  € 

where  Uk  snd  are  K-dimenslonal  spaces  of  C^IS^)  functions.  We  then 
multiply  (1.1)  by  V,  replace  u  by  U,  integrate  over  S^f  and  integrate 
the  time  derivative  and  diffusive  terms  by  parts  to  obtain  the  follow¬ 
ing  marching  problem  for  determining  U(x,t)  in  successive  strips 
S|^ ,  n  ™  0, 1 , . . .  • 

(2.2a)  U(x,0)  =  Puo(x),  a  <  X  <  b,  n  =  0, 


(2.2b) 


1  ^ 

F(V,U)  :=  /  J  {-V'^U  +  V^f (x,t,U,Ux)  +  v'rD(x,t,U)U  }dxdt 


^  ^n+1  ^n+1  ^ 

+  /  v'l^Udx  I  -  /  vTD(x,t,U)U  dt  I  =  0, 


Vv  e  1/k,  (x,t)  e  Sn,  n  >  0 


Here,  P  Is  an  interpolation  operator  on  the  space  and  U  must  also 
satisfy  any  essential  (Dirlchlet)  boundary  condition  in  (1.2b,c). 
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In  order  to  select  finite 


n 


into  N 

trapezoids  T^,  i  = 

,  n 

f  "  > 

vertices  ( 

c 

1 

(’'i '  ^n )' 

clement  b-ises  for  U^^  and  l/|^  we  partition 
n 

where  is  tJie  trapezoid  with 
n+1  ntl 

{’‘i-l'tn+1  j,  (Xi  ,  tn+i)  (cf.  Figure  1) 


Figure  1.  Space-time  discretization  for  time  ^■'tep  tn  <  t  < 


We  write  U  €  as 

K 

(2.3)  U(x,t)  =  Z  ci.(t)4.i(x,t) 

i  =0 

n  n 

where  each  (^j^(x,t)  is  selected  to  be  nonzero  only  on  (J  “^1  • 
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n 

Specifically,  we  map  each  Tj^  in  tl^e  (x,t) -plane  into  the  rectangle 
(2.4)  R  =  {(C,T)1-1  <5<1,  0<T<1} 


in  the  (C,T)-plane  and,  at  present,  we  choose  (^j(x,t),  j  =  0,1,..., K 
to  be  either  piecewise  C®  linear  or  a  piecewise  Hermite  cubic 

n 

polynomials  in  C  on  We  also  select  ^j(x,t),  i  =  0,1,..., K  as  a 

basis  for  l/^;  thus,  the  dimension  of  and  is  either  N  or  2N  for 
linear  or  cubic  approximations,  respectively. 

The  integrals  in  Eq.  (2.2b)  are  transformed  element-by-element  Into 
integrals  over  R  and  are  evaluated  numerically.  We  use  the  Trapezoid¬ 
al  rule  to  evaluate  the  T  integrals  and  a  three-point  Gauss -Legendre 
rule  to  evaluate  the  ^  integrals.  The  resulting  system  of  nonlinear 
algebraic  equations  is  solved  by  Newton’s  method,  with  users  supplying 
formulas  for  the  Jacobians  f^,  fy  ,  and  Dy.  Additional  details  on  our 

X 

finite  element  discretization  may  be  found  in  Davis  and  Flaherty  [7]. 


3.  Adaptive  Mesh  Selection  Strategies.  In  this  section  we  discuss 
several  algorithms  for  moving  the  mesh  so  that  the  spatial  discretiza¬ 
tion  error  in  L2  is  approximately  minimized  at  t  =  It  is  known 

(cf.,  e.g.,  (171)  that  the  spatial  error  in  finite  element-Galerkin 
methods  for  problems  like  (1.1.2)  satisfies  an  estimate  of  the  form 


(3.1)  lu-UI  <  Clu-Pul  , 

where  Pu  f  Uk  interpolates  the  solution.  If  we  assume  that  u(x,t)  € 
C^ta,bl,  then  Pereyra  and  Sewell  (16]  show  that  the  Interpolation 
error  in  (3.1)  can  be  asymptotically  minimized  in  L2  for  piecewise 
polynomial  interpolants  of  degree  k  -  1  by  selecting  the  Interpolation 
points  Xi(t),  i  =  0,1, .../N,  at  time  t  such  that 

k 

(3.2)  hi(t)  g(<:t,t)  =  E(t),  i  =  1,2, ..,,N  . 

Here, 

1  /2 

(3.3a,b)  h^lt)  =  Xi(t)  -X£_i(t),  g(x,t)  =  (x,  t)  [u(’t)  (x,  t)  ]  }  , 

u^*^)  is  the  kth  derivative  of  u  with  respect  to  x,  €  (xj^_i,  xj^),  E(t) 

n 

is  an  undetermined  function  of  t,  and  Xi(t)  is  the  line  joining  Xj^  and 
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n+1 

Xj  (cf.  Figure  1). 

The  result  (3.2)  states  that  the  interpolation  error  is  asymptoti¬ 
cally  minimized  when  the  mesh  is  moved  so  as  to  equldis tribute  the 
local  spatial  discretization  error.  Many  computer  codes  for  two-point 
boundary  value  problems  use  equidis tribution  algorithms  to  adapt  their 
computational  mesh  (cf.  Childs  et  al;  14)).  Additional  success  has 
been  reported  in  using  equidis  tribution  algorithms  for  variable  )tnot 
spline  interpolation  (cf.,  e.g.,  de  Boor  (8)). 

In  [7],  Davis  and  Flaherty  solved  (3.2)  by  an  iterative  technique. 
Herein  we  discuss  an  alternate  procedure  which  is  more  restrictive, 
but  has  a  simpler  structure.  To  begin,  we  follow  de  Boor  [8)  and 
ta)ce  the  Icth  root  of  (3.2)  and  write  it  in  the  asymptotically  equiva¬ 
lent  form 

1/k 

(3.4)  /  g(x,t)  dx  =c(t), 

Xi-1 

where  c(t)^  *■  E{t).  We  let 
X 

(3.5)  T(x,t)  =  /  g(s,t)Vk  ds  . 

a 

Then 


(3.6)  c(t)  =  (1/N)T(b,t) 

and  the  equidis tributtng  mesh  x^,  i  =  0,1,..., N,  is  determined  as  the 
solution  of  the  nonlinear  system 

(3.7)  T(xi,t)  =  ic(t)  ,  i  =  0,1, ...,N  . 

Of  course,  u^^J  is  unitnown  and  it  must  be  approximated  by  U.  To 
this  end,  suppose  that  we  have  computed  a  finite  element  solution 

n 

U(x,t,^)  at  time  t^  and  on  the  mesh  x^,  t  =  0,1,..., N.  We  differen¬ 
tiate  U(x,tn)  once  for  piecewise  linear  approximations  or  thrice  for 
Hermite  cubic  approximations  and  find  piecewise  constant  approxima¬ 
tions  for  0'(x,t„)  or  0''’(x,t^),  respectively.  We  then  use  finite 
difference  approximations  of  these  derivatives  to  compute  values  of 
u(^J(xi,tn)  and  g(xj^t„)  (cf.  (3.3b)),  for  I  =  0,1, ...,N  and  Jc  =  2  or 
4.  We  experimented  with  three,  four,  and  five  point  difference  formu¬ 
las  and  found  that  the  five  point  formulas  gave  marginally  better 
results,  so  we  are  using  them  in  the  current  version  of  our  code. 
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We  further  assume  that  g(x,t)V**,  k  =  2  or  4,  is  a  piecewise  linear 

n 

function  of  x  with  respect  to  the  mesh  x^,  i  =  0,1,..., N,  and 
integrate  it  to  find  a  piecewise  parabolic  approximation  to  T(x,t^) 
f  rom  (3.5). 

Finally,  we  find  cCt^)  using  (3.6)  and  determine  an  approximate  equi- 

..n 

distributed  mesh  Xj^,  i  =  0,  1,...,N,  at  time  t^  by  solving  (3.7)  using 
the  quadratic  formula. 

The  entire  process  can  potentially  be  iterated  to  find  a  better 
mesh;  however,  given  all  of  the  approximations  made  in  evaluating 

,  t,j ),  we  have  only  tried  this  at  t  =  0,  where  the  initial 
function  U()(x)  is  known  to  arbitrary  precision. 

The  equidlstrlbution  algorithm  has  a  nonunique  solution  whenever 
g(x,t)  =  0;  therefore,  we  may  expect  numerical  difficulties  whenever 
g(x,t)  is  small  on  any  subinterval.  We  combat  this  problem  by 
imposing  a  lower  bound  on  g,  i.e.,  we  replace  g(x, t)  by 

(3.8)  g(x,t)  ;=  g(x,t)  +  n  , 

in  Eqs.  (3.4)  and  (3.5).  Here,  n  is  a  small  empirically  determined 
quantity  that  is  discussed  further  in  Davis  and  Flaherty  [7).  Among 
other  things,  n  Insures  that  the  solution  of  (3.7)  is  a  uniform  mesh 
whenever  g(x,t)  is  small  everywhere  on  {a,b]. 

Our  discussion,  thus  far,  has  concerned  the  computation  of  an  equi- 
distributing  mesh  at  time  level  t^  where  a  solution  U(x,tn)  has 
already  been  computed.  To  obtain  an  estimate  for  an  optimal  mesh  at 
time  tn+1  prior  to  computing  the  solution  there,  we  extrapolate  the 
optimal  grids  from  previous  time  levels.  At  the  present  time,  we  are 

n  +  1  ,.n 

using  zero  order  extrapolation,  i.e.,  x^  =  xj^  ,  i  =  0, '1,...,N; 
however,  we  are  experimenting  with  several  different  extrapolation 
strategies,  some  of  which  are  discussed  in  Section  3.1. 


3.1.  Adaptive  Strategies  for  Mesh  Velocities.  The  zero  order 
extrapolation  strategy  that  was  just  discussed  was  applied  to  several 
examples  (cf.  Section  4)  and,  despite  its  simplicity,  it  has  worked 
quite  well,  even  on  problems  with  rapidly  moving  wave  fronts.  Never¬ 
theless,  we  can  expect  that  there  will  be  some  problems  where  it  will 
fall  to  produce  an  acceptable  mesh.  Most  of  our  attempts  to  use  higher 
order  extrapolation  produced  grids  that  oscillated  wildly  from  time 
step-to-timc  step,  even  when  the  solution  changed  quite  little.  In  an 
attempt  to  understand  and  remedy  this  phenomenon  while  simultaneously 
developing  a  more  dynamic  adaptive  mesh  strategy,  we  differentiated  Eq. 
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(3.4)  with  respect  to  time  and  obtained  the  followinq  system  for  the 
mesh  velocities: 

Xi 

(3.9)  Xi  q(Xi,t)  -  Xi_i  g(Xi_i,t)  +  /  gt(x,t)dx  =  c  , 

where  (  )  ;=  3(  )/3t.  This  system  offers  several  advantages  when  it 
is  used  in  conjunction  with  our  finite  element  scheme.  For  example, 
we  can  estimate  U^(x,t^)  while  assembling  the  finite  element  equa- 

n 

tions.  Then,  assuming  that  x^,  i  =  0,1, ...,N,  is  an  equidlstrlbuting 
mesh,  we  can  approximate  gt^^i't^)*  ^  =  0,1,..., N,  using  the  same 
finite  difference  scheme  that  was  used  to  find  gCx^.tj,),  i  =  0,1,..., N. 

Having  done  this,  c  is  determined  as 

(3.10)  c  =  (1/K)  /  gt(x,tn)dx  , 

a 

.n 

and  the  mesh  velocities  x^  ,  i  =  0,1, ...,N,  follow  readily  from  (3.9) 

n 

by  a  procedure  similar  to  the  one  that  we  used  to  find  x^  ,  i  = 

0,1,..., N.  We  can  then  integrate  the  mesh  velocities  using  an  explicit 
method  for  ordinary  differential  equations  and  obtain  an  approximation 

n+1 

for  an  equidlstrlbuting  mesh  ,  i  =  0, 1 , . . . ,N. 

However,  since  our  experiments  with  higher  order  (multi-level)  mesh 
extrapolation  produced  unstable  results  and  since  this  type  of  extra¬ 
polation  can  be  viewed  as  a  consistent  numerical  approximation  to  the 
differential  system  (3.9),  we  examine  the  stability  of  (3.9)  before 
proceeding  further.  Our  analysis  is  quite  general  and  is  not  limited 
to  either  the  specific  form  of  g(x,t)  that  is  given  in  (3.3b)  or  to 
piecewise  linear  mesh  trajectories. 

We  assume  that  xj^(t),  i  =  0,1, ...,N,  is  an  equidistributing  mesh 
that  exactly  satisfies  (3.4)  and  (3.9)  and  introduce  a  small  pertur¬ 
bation  6x^(0),  i  s:  0,1,..., M,  at  t  =  0  that  satisfies 

Xi+6xi 

(3.11)  /  g(x,0)dx  =  c(0)  +  6ci(0),  i  =  0,1, ...,N  . 

Xi-i+6xi_i 

Since  xq  and  xjj  are  fixed,  the  perturbations  must  also  satisfy 

N  N 

(3.12)  «xo(t)  =  6xN(t)  =  0  ,  E  «xi(t)  =0,  E  6ci(0)  =  0  . 

1=0  i=0 


e 


(3.16) 
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The  system  (3.9)  is  stable  to  linear  perturbations  when  6x(t) 
decays  and  this  occurs  when  •L~^{t)L(0>*  <  1/  for  some  nvatrix  norm. 
Unfortunately,  the  choice  of  g(x,t)  given  by  Eq.  (3.3b),  and  other 
reasonable  choices,  are  likely  to  be  decreasing  functions  of  time  for 
dissipative  parabolic  partial  differential  equations  and  this  will 
almost  certainly  yield  a  value  for  IL~^{t)L(0)l  that  is  larger  than 
unity. 

Local  instabilities  can  also  occur  when  the  mesh  is  moved  so  that 
one  of  the  three  ratios  Involving  g  in  Eq.  (3.16)  exceeds  unity.  How¬ 
ever,  since  these  instabilities  are  local,  they  may  either  grow  or 
decay  as  time  progresses . 

The  following  two  examples  Illustrate  some  of  the  instabilities  that 
can  occur  in  Eq.  (3.9). 

Example  3.1.  Consider  the  constant  coefficient  heat  conduction 
problem 

(3.17a)  Ut  =  Ujjjj  ,0<x<1,t>0, 

(3.17b,c)  u(x,0)  =  simix  ,  u(0,t)  =  u(1,t)  =  0, 


which  has  the  exact  solution 

(3.17d)  u(x,t)  =  e'*’’  ^  sinirx  . 

Since  the  solution  of  this  problem  Is  separable,  the  optimal  strategy 
is  to  generate  an  equidistributed  mesh  at  time  t  =  0  and  use  it  for  all 

2 

time.  When  g(x,  t)  is  given  by  (3.3b)  we  find  that  lL~^(t)L(0)l  -e"’^ 
thus,  we  expect  the  solution  of  (3.9)  to  t)e  unstable.  In  Figure  2,  we 
display  the  meshes  produced  by  both  (3.4)  and  (3.9)  for  g(x,t)  = 
|uyx(x,t)|,  and  the  instability  produced  by  using  (3.9)  is  clearly 
visible.  We  note  that  exact  values  of  u^x  were  used  in  obtaining 
Figure  2  and  the  only  errors  that  were  introduced  in  the  computation 
were  due  to  trapezoidal  rule  integration  and  a  perturbation  of  0.01 
in  the  initial  mesh  for  Eq.  (3.9). 

Example  3.2.  We  consider  a  constant  coefficient  heat  conduction 
problem  that  was  studied  in  Davis  and  Flaherty  (7),  i.e., 

(3.18)  Ut  =  Ujjx  +  f(x,t)  ,0<x<1,t>0. 

The  initial  conditions,  Dirichlet  boundary  conditions,  and  source  f 
are  chosen  so  that  the  exact  solution  is 

(3.19)  u(x,t)  =  tanh(ri(x-1)  +  r2t)  . 
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S4 _ H _ U - 1L_*-U U ILhJI — u — »u - 1 

iun  O)  0.40  X  160  0.80  UO 


Figure  2.  Meshes  that  were  produced  by  solving  Eq.  (3.4)  (solid  line) 
and  integrating  Eq.  (3.9)  (dashed  line)  for  Example  3.1. 


The  solution  (3.19)  is  a  wave  that  travels  in  the  negative  x  direction 
when  ri  and  V2  are  positive.  The  values  of  r^  and  r2  determine  the 
steepness  of  the  wave  and  its  propagation  speed.  The  meshes  produced 
by  both  (3.4)  and  (3.9)  for  r^  =  r2  “  5  and  g(x,t)  = 

shown  in  Figure  3.  The  solution  of  Eq.  (3.9)  Is  Initially  stable,  but 
as  time  progresses  and  g(x,t)  decays  it  becomes  unstable. 
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Figure  3.  Meshes  that  were  produced  by  solving  Eq.  (3.4)  (solid  line) 
and  integrating  Eq.  (3.9)  (dashed  line)  for  Example  3.2. 


Petzold  (161  suggested  that  the  following  linear  combination  of 
Eqs.  (3.4)  and  (3.9)  might  yield  stable  meshes  with  some  improved 
dynamic  behavior: 


Xig(xi,t)  -  Xi_ig(xi_i,t)  +  }  gt(x,t)dx  +  X  j  g(x,t)dx 

xi-1  Xi_, 


(3.20) 


c  +  Xc 
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Here  A  >  0  is  a  parameter  to  be  determined  so  that  (3.20)  is  stable. 

We  have  not  done  enough  analysis  or  experimentation  to  form  any  firm 
conclusions  ;  however.  In  Figure  4  we  compare  the  solution  of  Eg.  (3,4) 
with  Eq.  (3.20)  for  Example  3.2.  Equation  (3.20)  was  solved  by  the 
explicit  Euler  method  using  A  =  1/At,  uniform  time  steps  of  At  =  1/20, 
exact  values  of  g  and  g^.  and  trapezoidal  rule  Integration.  It  appears 
tliat  the  mesh  produced  by  Eq.  (3.20)  has  some  local  instabilities,  but 
is  globally  stable  for  this  example.  Hence,  this  method  has  promise, 
but  much  more  testing  is  needed  to  determine  its  behavior,  especially 
when  the  approximate  solution  U(x,t)  Is  used  to  calculate  g  and  g^. 


Figure  4.  Meshes  that  were  produced  by  solving  Eq.  (3.4)  (solid  line) 
and  integrating  Eq.  (3.20)  (dashed  line)  for  Example  3.2. 


4.  Examples .  In  this  section  we  examine  the  performance  of  our 
adaptive  finite  element  method  on  four  examples.  In  all  of  these  we 
used  the  equidistribution  algorithm  as  discussed  in  Section  3  with 
zero  order  extrapolation. 
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Example  4.1.  We  consltler  Example  3.2  (cf.  Bqs.  (3.18)  and  (3.19)) 
with  ri  =  r2  =  100.  We  compare  the  polntwise  errors  vs.  x  at  t  =  1  fi 

computations  performed  with  uniform  an«l  adaptive  meshes  of  20  finite 

elements  using  piecewise  linear  (cf.  Figure  5a)  and  cubic  approxima¬ 
tions  (cf.  Figure  5b).  The  results  were  obtained  using  uniform  time 
steps  of  At  =  0.01  and  0.0025,  respectively,  for  the  linear  and  cubic 
approximations.  We  see  that  the  computations  using  the  fixed  uniform 
grids  have  large  errors  at  the  wave  front,  which  is  near  x  =  0  at  t  = 

1.  The  adaptive  mesh  algorithm,  on  the  other  hand,  concentrates 

elements  In  the  wave  front  (cf.  Figure  3)  and  distributes  the  local 
error  more  evenly  over  the  domain. 


Figure  5.  Local  error  at  t  =  1.0  for  Example  4.1  with  20  elements. 

The  solid  curve  was  computed  on  a  fixed  uniform  mesh  and  the  dashed 
curve  on  an  adaptive  mesh.  The  solution  in  Figure  5a  (left)  used 
piecewise  linear  approximations  with  uniform  time  steps  of  At  =  0.01 
and  that  in  Figure  5b  (right)  used  piecewise  cubic  approximations  with 
uniform  time  steps  of  At  =  0.0025. 


Example  4.2.  We  consider  the  following  problem  for  Burgers'  equa¬ 
tion: 

u^  +  uUjj  =  »0<x<1,t>0 

(4.1) 

u(x,0)  =  sinrrx  ,  u(0,t)  =  u(1,t)  =  0  , 


and  c  =  5x10"^.  The  solution  of  this  problem  is  a  pulse  that  steepens 
as  it  travels  to  the  right  until  it  forms  a  shock  layer  at  x  «  1. 

After  a  time  of  0(1/e)  the  pulse  dissipates  and  the  solution  decays  to 
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zero.  We  compare  solutions  using  piecewise  Hermite  cubic  approxima¬ 
tions  on  a  uniform  and  an  adaptive  mesh  of  10  finite  elements  in 
Figures  6a  and  6b,  respectively.  Both  calculations  were  performed  with 
uniform  time  steps  of  At  =  0.1.  It  is  well  known  that  finite  differ¬ 
ence  and  piecewise  linear  finite  element  solutions  of  this  problem 
exhibit  spurious  mesh  oscillations  unless  the  mesh  width  is  0(e) 
within  the  shock  layer.  However,  the  solution  using  piecewise  cubic 
approximations  on  a  uniform  mesh,  that  is  shown  in  Figure  6a  at  t  = 

0.6,  is  pointwise  very  accurate  and  the  large  errors  appear  in  the 
slope  of  the  computed  solution.  These  errors  largely  disappear  when 
the  mesh  adapts  with  the  solution  and  is  appropriately  concentrated  in 
the  shock  layer  (cf.  Figure  6b). 


Figure  6.  Solution  of  Example  4.2  using  cubic  approximations  with  10 
elements  and  uniform  time  steps  of  At  =  0.01.  The  solution  in  Figure 
6a  (left)  was  obtained  using  a  fixed  uniform  mesh  while  that  in  Figure 
6b  (right)  used  an  adaptive  mesh. 


Example  4.3.  We  are  currently  investigating  the  following  focusing 
problem  for  the  nonlinear  Schrodinger  equation  in  cylindrical  coordi¬ 
nates: 

(4.2)  ut  *  (i/r)(ruj.)j.  +  i|u|2u  ,  r  >  0  ,  t  >  0  ,  u(r,0)  =  Ae“®*^  , 

where  i  =  and  u(r,t)  is  a  complex  valued  function. 
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It  is  known  (cf.  (121)  that  the  solution  of  (4.2)  “self-focuses", 
i.e.,  it  develops  infinite  values  of  u  in  a  finite  time,  when  the 
initial  conditions  are  "strong  enough".  Problems  of  this  type  occur 
in  laser  optics,  and  we  are  trying  to  determine  the  conditions  which 
cause  the  solution  to  "blow  up"  and  also  the  local  character  of  the 
solution  just  prior  to  blow  up.  This  problem  is  still  under  investi¬ 
gation,  in  collaboration  with  A,  Newell  15],  and  herein  we  only 
present  some  preliminary  results  which,  we  feel.  Illustrate  the  need 
for  adaptive  mesh  strategies  on  difficult  nonlinear  partial  differen¬ 
tial  equations. 

In  Figures  7a  and  7b  we  illustrate  ju(r,t)|  for  two  sets  of  initial 
conditions  having  a  =  25  and  A  =  14  and  A  =  15,  respectively.  The 
solution  in  Figure  7a  does  not  focus  while  (we  speculate  that)  the 
solution  in  Figure  7b  focuses.  The  dramatic  growth  in  the  magnitude 
of  the  solution  that  is  shown  in  Figure  7b  is  also  accompanied  by 
rapid  changes  in  phase  as  focusing  nears. 


Example  4.4.  We  are  studying  elastic-plastic  impact  problems  for 
cylindrical  rods  using  the  long  rod  model  of  T.  W.  Wright  [18]: 


(4.3)  V  =  Wt  ,  e  =  Wjj  ,  P  *  ut  »  q  =  Ux  , 


(4.4)  et  =  Vx  ,  qt  =  Px  ' 


(4.5)  v^  *  Sx  ,  p^  —  2(Qx— P)  . 


Here,  u  and  w  are  dimensionless  radial  and  axial  displacement  compon¬ 
ents,  p  and  V  are  radial  and  axial  velocity  components,  e  is  the  axial 
strain,  q  is  the  shear  strain,  and  S,  P,  and  Q  are  axial,  radial,  and 
shear  forces,  respectively.  Equations  (4.3)  define  the  strain  and 
velocity  variables  in  terms  of  the  displacements,  Bqs.  (4.4)  are  com¬ 
patibility  relations,  and  Eqs.  (4.5)  are  the  equations  of  motion.  We 
also  need  appropriate  constitutive  laws  that  relate  S,  P,  and  Q  to  e, 
u,  and  q,  and  herein  we  simply  use  the  linear  Hooke’s  laws 


(4.6)  S  =  e  +  u  ,  P  =  2  (ve+u)  ,  Q  =  -  q  , 

1-v  1-v  4(  1-v) 


Figure  7«  Magnitude  |u|  of  the  solution  of  the  Schr^inger  equation 
(4.2)  with  a  =  25  and  A  »  14  (Figure  7a,  left)  and  A  ^  15  (Figure  7b, 
right) . 


where  v  is  Poisson's  ratio.  If  v  =  0  we  have  a  one -dimensional  theory 
and  Eqs.  (4.3)-(4.6)  reduce  to  a  simple  wave  equation.  However,  if 
V  ^  0  Bqs.  (4.3)-(4.6)  give  a  two-dimensional  theory  that  is  valid 
when  the  length  L  of  the  cylindrical  rod  is  large  compared  to  unity. 

Lilce  the  previous  example,  this  is  work  in  progress  in  collabora¬ 
tion  with  J.  J.  Wu  and  T,  W.  Wright,  so  we  will  only  report  some  pre¬ 
liminary  results  that  illustrate  the  differences  between  the  one-  and 
two-dimensional  theories.  We  use  homogeneous  initial  conditions  and 
the  boundary  conditions 


(4.7)  v(0,t)  -  0.1  ,  Q(0,t)  =  0  ,  S(5,t)  »  0  ,  Q(5,t)  »  0 
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These  conditions  correspond  to  a  cylinder  of  length  L  =  5  that  is  hit 
at  t  =  0  by  a  rigid  wall  that  is  moving  with  a  velocity  of  ten  percent 
of  the  longitudinal  wave  speed  of  the  bar.  Since  the  initial  condi¬ 
tions  are  homogeneous,  bur  adaptive  algorithms  have  no  choice  but 
to  select  a  uniform  initial  mesh.  However,  the  velocities  and  strains 
jump  at  X  =  0,  t  =  0  and  the  mesh  should  be  concentrated  in  the  vicin¬ 
ity  of  this  point.  There  are  several  possibilities  for  overcoming 
this  difficulty.  For  example,  we  could  either  Input  a  graded  Initial 
mesh  or  we  could  use  nonhomogeneous  initial  conditions  and  start  the 
problem  at  a  small  value  of  t  >  0.  We  tried  both  alternatives,  and 
they  gave  very  similar  results  for  t  >  0.  We  also  added  artificial 
viscosity  terms  ce^x'  ^'^xx'  ^Pxx»  respectively,  to  the  right 

sides  of  Eqs.  (4.4)  and  (4.5). 

We  present  results  for  the  axial  stress  S  vs.  x  for  t  =  0,1,2,. ..,10 
In  Figures  8  and  9  for  v  =  0  and  v  =  0.3,  respectively.  In  each  case 
50  piecewise  linear  elements  were  used  with  an  artificial  viscosity 
coefficient  €  =  0.005. 

It  Is  known,  that  one -dimensional  elastic  waves  are  non-dlspersive 
whereas  two-dimensional  waves  are  dispersive.  The  dispersive  nature  of 
the  two-dimensional  solution  is  clearly  visible  in  Figure  9.  Unfortu¬ 
nately,  the  one -dimensional  solution  is  much  more  dispersive  than  It 
should  be.  We  conjecture  that  the  excessive  dispersion  is  due  to  our 
form  of  the  artificial  viscosity,  and  we  are  experimenting  with 
different  models. 

5.  Discussion.  The  computations  of  Section  4  show  that  it  is 
possible  to  construct  an  accurate  and  stable  adaptive  grid  finite 
element  procedure  for  nonlinear  systems  of  partial  differential  equa¬ 
tions  that  offers  several  advantages  over  fixed  grid  methods.  However, 
we  have  much  more  to  do  in  order  to  achieve  our  goal  of  developing 
reliable  and  robust  general  purpose  software  for  partial  differential 
equations . 

In  this  paper  we  concentrated  on  improving  the  mesh  moving  capabil¬ 
ities  of  our  code.  We  have  always  found  equldlstrlbution  with  zero 
order  extrapolation  to  be  unsatisfying;  however,  our  attempts  to  use 
higher  order  extrapolation  always  ended  in  failure.  The  theoretical 
results  of  Section  3  explain  why  this  is  so,  and  Indicate  a  possible 
remedy  that  we  hope  to  explore  further. 

We  feel  that  we  have  reached  the  limit  in  terms  of  what  can  be 
achieved  with  a  fixed  number  of  elements  per  time  step.  There  are 
basically  two  possible  ways  of  extending  our  methods  to  include  the 
ability  of  adding  and  deleting  elements  as  the  integration  progresses 
in  time.  We  can  envision  a  method  of  lines  approach,  in  the  spirit  of 
Bieterman  and  Babuska  [2,3];  however,  with  "lines"  that  adaptively 
equidlstribute  the  local  spatial  component  of  the  error.  Elements  can 
be  added  or  deleted  as  the  Integration  progresses  and  the  power  of  the 
ordinary  differential  equations  codes  can  still  be  used. 
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A  second  possibility  is  to  locally  refine,  using  space-tlme-trape- 
zoidal  elements  in  regions  having  large  local  error.  This  Is  In  the 
spirit  of  the  adaptive  finite  difference  methods  of  Berger  Ml;  how¬ 
ever,  the  computational  cells  would  be  trapezoidal  rather  than 
rectangular. 

We  are  exploring  both  approaches,  but  feel  that  the  latter  offers 
the  most  promise. 
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